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Summary  i 

* i 

Two  principal  and  related  topics  are  considered:  (1)  adaptive  mesh 
refinement  for  finite  element  computations,  and  (2)  mesh  refinement  specifi- 
cally for  nonlinear  flow  problems.  In  the  first  instance  the  residual  and 
&sociated  trace  theorems  for  variational  problems  are  introduced  to  relate 
the  solution  error  to  a computable  residual.  This  provides  a theoretical 
basis  for  a mesh  refinement  strategy. 

A corresponding  adaptive  refinement  procedure  that  automatically  and 
selectively  refines  the  mesh  is  formulated  and  implemented  for  a class  of 
nonlinear  transport  problems  in  chemical  engineering.  In  particular,  a non- 
linear problem  with  a boundary-layer  solution  is  investigated.  The  strategy 
of  interweaving  Newton  solution  and  mesh  refinement  proves  particularly  ef- 
ficient. 

Two-dimensional  compressible  and  transonic  flows  are  next  examined. 

Mesh  refinement  of  subregions  of  the  flow  field  is  applied  to  yield  high 
solution  near  a singularity  for  both  the  linear  and  nonlinear  flows.  Re- 
finement and  Newton  iteration  are  combined,  together  with  Mach  number  para- 
meterization, to  determine  an  efficient  and  accurate  solution  algorithm. 

Similar  points  for  nonlienar  viscous  problems  are  also  reviewed. 

_ 

Invited  lecture  presented  at  FENOMECH  '78  (Finite  Elements  and  Nonlinear 
Mechanics),  Stuttgart,  and  to  appear  in  the  conference, proceedings. 
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1.  Introduction 

In  approximating  the  solution  of  boundary-value  problems  by  finite 
element  methods,  we  are  continually  confronted  with  the  problem  of  deciding 
on  a mesh,  assessing  the  quality  of  the  solution  obtained  with  that  mesh, 
and  perhaps  iteratively  adjusting  and  enriching  the  mesh.  As  the  solution 
and  its  derivatives  may  vary  markedly  over  the  domain,  nonuniform  meshes 
are  warranted  in  the  interests  of  computational  efficiency.  This  is  es- 
pecially true  of  nonlinear  problems  where  solutions  exhibit  layer  or  singu- 
lar behavior. 

Of  course,  this  refinement  problem  applies  quite  generally  to  discrete 
methods  for  boundary-value  problems.  In  the  present  article  we  shall  confine 
the  choice  of  applications  to  flow  problems.  Fluid  mechanics  has  long  been 
a rich  source  of  problems  in  applied  mathematics.  Many  classical  methods 
of  analysis,  such  as  analytic  function  theory,  have  their  origins  and  motiva- 
tion in  fluid  mechanics.  Problems  with  singularities,  interior  layers,  boun- 
dary layers  and  nonlinearities,  where  meshes  are  important,  are  commonplace. 
This  is  true  of  the  entire  gamut  of  viscous,  compressible,  rotational  and 
other  flows  and  also  of  transport  problems  which  are  included  here. 

The  focus  of  this  article  is  the  underlying  mathematical  analysis,  fi- 
nite element  technique,  algorithms  and  implementation  of  mesh  refinement  in 
this  context.  Special  points  of  interest  concerning  refinement  in  layers 
and  near  singularities  , and  solution  strategies,  are  developed  for  representa- 
tive steady  flow  problems.  Throughout,  important  points  are  demonstrated  using 

1 


— ^ — .... 


2 


select  problems  drawn  primarily  from  ongoing  research  studies  of  this  topic. 
Similar  investigations  have  also  been  carried  out  for  nonlinear  problems  in 
solid  mechanics  and  work  in  both  areas  is  continuinp.  Where  appropriate, 
other  research  activities  related  to  this  topic  are  also  indicated. 

2.  Refinement  and  Residuals 

In  many  engineering  applications  one  is  able  to  predict  those  subregions 
where  the  solution  is  changing  markedly  and  use  a finer  mesh  there.  For  in- 
stance, the  form  of  the  loading  or  boundary  data  and  shape  of  the  domain 
may  indicate  where  the  mesh  needs  to  be  fine.  A preliminary  computation  on 
a coarse  uniform  mesh  may  also  be  useful.  Yet  these  actions  only  alleviate 
the  situation  and  do  not  resolve  the  main  problem  - the  degree  of  local  sub- 
division is  still  unknown. 

Ideally,  we  require  an  algorithm  which  selectively  refines  the  mesh  in 
appropriate  subregions  of  the  domain.  This  implies  that  we  have  some  cri- 
terion for  automatically  assessing  the  quality  of  the  approximate  solution 
on  intermediate  finite  element  meshes.  It  is  easy  to  devise  ad  hoc  princi- 
ples for  refining  the  mesh.  Monitoring  the  changes  in  local  solution  beha- 
vior between  successive  meshes  and  extrapolating  to  the  next  mesh  is  perhaps 
the  simplest  strategy  and  will  often  work  reasonably  well.  We  also  note 
that  one  can  readily  concoct  pathological  examples  such  as  highly  oscillatory 
functions  that  will  defy  any  refinement  procedure  (Try  u = — sin  ex, 

0 < e « 1). 

A criterion  which  has  a more  rigorous  foundation  in  the  approximate 


method  itself  is  more  desirable.  Finite  element  methods  are  often  based  on 
energy  principles:  the  energy  functional  is  minimized  for  the  exact  solution, 
and  the  approximate  solution  minimizes  the  same  functional  over  a finite 
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dimensional  subspace  S(M),  the  class  of  finite  element  aooroximations  on  a given 
mesh  M . Refining  the  mesh  produces  a broader  class  of  functions  and  a 
lower  value  to  the  approximate  global  energy.  This  suggests  that  computation 
of  the  energy  for  successive  meshes  will  indicate  the  global  quality  of  the 
solution  and  meshes.  Perhaps  one  can  extend  this  argument  to  the  local  ac- 
curacy of  the  energy  contribution  from  individual  elements  or  groups  of  ele- 
ments as  they  are  refined. 

What  may  we  infer  regarding  those  problems  where  a classical  energy 
principle  is  not  available?  The  weak  variational  statement  and  generalized 
solution  constitute  an  equivalent  integral  formulation  of  a boundary- value 
problem.  The  classical  minimization  principles  are  then  included  as  a spe- 
cial category.  In  turn,  we  recognize  that  the  variational  problem  can  be 
interpreted  as  a particular  form  of  the  method  of  weighted  residuals.  We 
now  appeal  to  this  characterization  to  establish  a suitable  refinement  cri- 
terion. For  a given  approximation,  the  residual  represents  the  amount  by 
which  the  differential  equation  fails  to  be  satisfied.  An  exact  classical 
solution  yields  residual  R = 0.  The  variational  problem  is  constructed  from 
a duality  pairing  between  the  residual  R(u)  and  test  functions  v . 

The  various  approximation  methods  require  that  the  residual  be  close  to 
zero,  in  some  sense.  That  is,  the  projection  of  the  residual  in  the  test 
space  is  zero.  For  example:  in  collocation  the  test  functions  are  delta 
functions  and  the  residual  is  forced  to  be  zero  at  specific  points;  in 
Galerkin  methods  the  projection  of  the  residual  in  the  approximation  space 
or  a subspace  is  zero  - in  an  appropriate  norm,  a point  we  shall  shortly 
develop-  This  leads  us  to  explore  the  use  of  the  residual  as  an  indication  of 
solution  accuracy  globally  in  the  domain  and  perhaps  even  locally.  Except  for 
very  simple  bases,  the  residual  will  usually  be  non-zero  and  well  defined  on  an 
element  and  the  norms  of  the  residual  easily 
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calculated  there.  For  example,  we  could  write  for  element  the  usual  l_2  norm 


Then  a refinement  scheme  might  be  introduced  that  refines  those  elements 

where  ||  R |L  per  unit  volume  is  large.  The  simplest  case,  a piecewise  linear 
“e 

element  draws  our  attention  to  a difficulty  here.  The  second  order  deri- 
vatives u"  will  not  contribute  to  the  residual  in  the  element  interior  but 
rather  as  boundary  terms  (jumps)  at  the  element  interfaces.  These  should  be 
included  in  the  residual  evaluation  of  eauation  (1). 

The  foregoing  discussion  motivates  a residual  analysis  to  provide  the 
desired  refinement  criterion.  In  principle,  we  seek  global  and  local  error 
bounds  for  the  solution  in  terms  of  computable  element  residuals.  We  shall 
develop  this  basic  analysis  for  the  standard  two-point  problem.  The  trace 
theorems  allow  us  to  employ  the  same  residual  analysis  in  higher  dimensions 
[11- 

Following  this  theoretical  analysis,  residuals  are  employed  in  adaptive 
refinement  and  related  procedures.  A particular  point  we  shall  examine  is  that 
of  interweaving  adaptive  mesh  refinement  with  nonlinear  solution  iteration.  Numeri- 
cal experiments  are  conducted  to  confirm  that  this  is  computationally  very  effective. 

The  particular  examples  described  concern  nonlinear  heat  and  mass  transport  with  a 
boundary  layer  solution  and  compressible  flow  problems. 

3.  One-Dimensional  Analysis 

The  one-dimensional  problem  has  a very  special  connectivity  - adjacent 
elements  meet  at  their  end  nodes.  Since  we  are  considering  mesh  refinement 
strategies  here,  this  implies  that  refinement  by  multiple  "splittings"  of 
elements  is  relatively  straightforward.  New  nodal  points  are  introduced  and 
the  continuity  and  differentiability  of  the  approximation  (C°,c\...  functions) 
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are  readily  maintained. 

We  first  summarize  some  fundamental  properties  of  the  residual  for 
linear  two-point  problems.  The  details  of  the  proofs  are  described  else- 
where [2-4].  Consider  the  standard  linear  problem, 

— (p(x)u* ) + q(x)u  = f(x)  in  fi=(0,l)  (2) 

with  u(0)  = u ( 1 ) = 0. 

A generalized  (weak)  solution  to  equation  (2)  can  be  determined  among 
admissible  functions  u(x) € Hg(ft)  by  requiring  that 


j pu'v'  dx  + j quv  dx  - J fv  dx  = 0 (3) 

0 0 0 

for  al 1 v(x)  6 Hg(n). 

Let  L denote  the  linear  differential  operator  in  equation  (2).  Then, 
for  an  approximate  solution  U(x)€  Hq(si)  the  residual  R(x)  = LU  - f is  de- 
fined in  a distributional  sense  by  equation  (3).  Since  Lu  = f and  LU  = R + f, 
then  L(U  - u)  = R or  Le  = R where  e = U - u is  the  solution  error.  Now  as  u 


J 


and  U are  in  Hg,  so  is  e , and  R is  in  the  dual  space  (Hg)‘  = H"^  . From 
Hilbert  Space  Theory,  the  induced  norn  in  H”^  is 


II  R II  , 


< R,w  > I 


w € 
w f 0 


Since  Le  = R with  Le  = -(pe1)'  + qe  , then  from  (3)  we  have  on  setting 
= e and  simplifying 


c Hi ,n  - CH  R 


where  C is  a constant. 

This  basic  inequality  forms  the  starting  point  of  our  residual  analysis. 
The  H”^-norm  is  not  readily  computed,  so  we  seek  to  replace  the  bound  by  an 
expression  involving  computable  element  residuals.  To  evaluate  ||  R||_-j  , 
we  return  to  the  definition  of  equation  (4)  and  examine  the  duality  pairing. 
On  a partition  of  n elements  of  C°-Lagrange  type. 


n (• 

<R,w>  = l (pu'w'  + quw  - fw)dx 

1=1  i 


(6) 


Integrate  by  parts  on  the  entire  domain  as  a union  of  elements  to  obtain 


+ f (pu'w)'  dx  | 

A J 


<R,w>  =J||  {- (pu  ‘ ) ' + qu  - f }w  dx  + ! (pu'w)'  dx  | (7) 

°i 

Taking  absolute  values  and  applying  the  Cauchy- Schwarz  inequality  to  each 
integral,  we  obtain 

| < R,w  > | < (||  R ||0f^  ||  w ||0jfi_  + hT^  |[pu'(xi+1)  - pu'(xi)]||  w \\un 

(8) 

Since  ||  w ||Q  < ||  w \\]  , 


< R,w>  | <_  {||  R ||0jfi  + h^|pu'(xi+1)  - pu  ’ (xi ) | } ||  w ||ljQ 


(9) 


From  the  Schwarz  and  Holder  inequalities  for  finite  sums. 


< R,w  > | < ||  R||0f^  + h“*s(|pu,(xi+1)  - pu,(xi)|)}2j^|| 


(10) 


We  define  the  specific  element  residual  as 
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R Mi  = Ml  R||0>n.  + h^(|Pu'(x1+1)  - pu,(xi)|)> 


(n) 


so  that 


| <R,w>  | <<I  «lRlll?>*ll«lli,n 

i=l 

Using  this  in  equations  (4.4), 

ii ril,  < H iii" hi -V - m*r 


(12) 


'n 


(13) 


and  inequality  (4.9)  is  replaced  by  the  global  residual  bound 


c Hi, a - c III R llln 


(14) 


The  element  contribution  to  |||R|||fi  may  be  computed,  within  a constant  C, 


from  equation  (13).  Since  C is  independent  of  h^  and  w , the  inequality 
(14)  is  a practical  bound  that  can  be  utilized  directly  in  our  refinement 
criterion.  As  |||R|||.  is  reduced  by  refinement,  |||R|||n  will  also  decrease 
as  will  the  global  error  ||  e 1 1 ^ ^ . 

To  obtain  an  effective  local  refinement  criterion  we  need  a local  bound 
on  an  element  of  the  same  form  as  inequality  (14).  Consider  a subsidiary 
problem  posed  on  an  element  and  with  nonhomogeneous  boundary  data  (approxima- 
tions to  the  solution  values)  at  the  element  end  nodes.  By  similar  arguments 
to  the  above  we  obtain  the  local  error  bound 


eHl,fie  idP 


'i.n„ 


(15) 


where  n arises  from  the  nonhomogeneous  data  and  satisfies  Ln  = 0 in  , 


n = e at  the  end  nodes  of  an  element.  Using  the  trace  theorems  and  maximum 


I 


V.  ~ j. 
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principles,  ||  n |U  n is  negligible  for  h = max  h.  sufficiently  small  and 
i ,»e  i i 

ftg  not  near  the  boundary  of  ft  . Hence,  inequality  (14)  provides  the  desired 
local  bound. 

4.  Algorithm 

This  residual  analysis  establishes  a theoretical  basis  for  developing 
an  adaptive  refinement  strategy  for  two-point  problems.  We  now  utilize  this 
residual  criterion  in  a statistical  procedure  in  which  subregions  are  located 
and  refined  automatically  as  part  of  the  algorithm.  The  statistical  approach 
is  quite  robust  as  we  demonstrate  in  numerical  experiments. 

Consider  an  initial  mesh  Mg(n).  A finite  element  approximation  U(x)  is 
computed  on  Mg(ft)  and  the  element  residuals  computed  according  to  equation 
(11).  The  mean  R^Mg)  and  standard  deviation  s(Mg)  are  calculated  for  the 
set  of  element  residuals  {|||R|||..}.  We  define  refinement  “intervals"  on  the 
distribution  of  element  residuals  by  considering  the  intervals  I-j  = (-<*>  , 

R + ks),  = (R  + ks,R  + 2ks ),...,  = (R  + (t-1  )ks,°°) , where  k represents 

the  fraction  of  a standard  deviation  in  each  interval.  For  example,  if  k = 1 
the  refinement  intervals  lie  l,2,3,...,t  standard  deviations  above  the  mean 
residual.  Now  if  |||R|||.  is  located  in  refinement  interval  j,  then  element  i 
is  refined  to  j elements.  An  estimate  of  the  solution  at  the  newly-intro- 
duced nodes  is  provided  by  local  interpolation  in  the  old  basis.  On  comple- 
tion of  the  refinement  of  Mg (ft)  to  M^(ft)  a solution  can  be  computed  on  M^(ft). 

The  residual  calculations  are  straightforward  and  readily  implemented 
in  a computer  program.  In  the  algorithm  above  we  begin  with  an  initial  coarse 
mesh  Mg(n)  and  the  algorithm  selectively  refines  the  mesh  in  a systematic 
manner  to  the  final  mesh  M^(fj)  with  solution  u(x;Mf(ft)). 

s not  important  that  the  criterion  be  precise  at  each  mesh  iteration. 
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Rather,  the  sequence  of  meshes  generated  need  only  approach  an  appropriate 
final  configuration.  The  means  {R}  of  successive  residual  distributions 
regress  towards  zero  and  a natural  gradation  of  the  mesh  is  inherent  in  the 
scheme.  These  points  are  demonstrated  in  the  example  below. 

5.  Example:  Interior-Layer  Problem 
Consider  the  two-point  problem 

i _ r>  i „ _ i _ 

-[{  1 + a(x  - x)  }u‘]'  * 2 [1  + a(  x - x)tan  u(x  - x)  + tan  ax) 

in  0 < x < 1 with  u(0)  = u(l)  = 0.  If  a is  large,  say  (0(10^)),  then  the 
solution 

u(x)  * (1  - x)[tan_1a(x  - x)  + tan  ux)  (17) 

has  an  interior  layer  in  the  neighborhood  of  x = x. 

Numerical  results  are  presented  here  for  a = 100  and  x = .36388.  In 
the  numerical  experiment  quartic  elements  are  used,  beginning  with  an 
initial  mesh  Mg(fl)  of  four  equal  elements.  After  8 adaptive  refinements,  the 
finite  element  mesh  consists  of  29  elements  with  a gradual  transition  from 
either  side  to  smaller  elements  within  the  layer  [7  in  (0,.25);  4 in  (.25, .325); 

9 in  (.325, .375);  4 in  (.375, .425);  5 in  (.425,1.0)].  At  this  mesh  level  the 
global  residual  is  7.479.  On  subsequent  refinements  many  elements  are  located 
in  the  interval  (0.25,0.5),  particularly  at  the  interior  layer.  On  the  "final" 

mesh  of  74  elements  (43  of  them  in  (.3125, .4375])  the  global  residual  is  .812. 

1 0 -3 

The  H and  H norms  of  the  global  solution  error  are  9.895  x 10  and  1.411  x 

-4 

10  , respectively. 


It  is  instructive  to  compare  global  residual  and  errors  for  the  entire 
refinement  history.  In  Figure  1 we  present  typical  results  as  a log-log  plot 
of  global  residual  |||R|||  versus  H°  and  norms  of  the  solution  error.  The 
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straight  line  corresponds  to  a linear  regression  of  these  results.  The  resi- 
dual and  error  decrease  monotonically  on  successive  refinements  and  the  ex- 
perimental values  are  consistent  with  earlier  theoretical  results.  Empiri- 
cally, we  have  for  this  example,  ||  e ||Q  = 4.703  * 10”6|||R  |||^ ,7^6  and 
||  e ||i  = 2.212  x 10"3 HI R H)1’517  . The  H1  error  norm  (slope  0.66)  decreases 
more  slowly  than  the  H°  error  norm  (slope  0.59).  Qualitatively,  this  re- 
sembles the  theoretical  result  for  uniform  mesh  refinement.  Data  points  for 
a uniform  mesh  refinement  are  also  plotted  and  lie  on  the  linear  regression 
curves  in  the  figure.  This  confirms  the  argument  that  the  ability  to  compute 
on  an  adaptive  non-uniform  grid  using  the  residual  criterion  essentially  en- 
ables one  to  progress  down  the  straight  lines  rapidly,  reducing  the  error  in 
the  solution  while  incurring  low  storage  and  computation  penalties. 

[Figure  1] 


6.  Nonlinear  Problems 

When  nonlinear  problems  are  considered,  automated  mesh  refinement  is 
even  more  appealing,  particularly  if  the  solutions  are  not  very  smooth  or  if 
layers  again  occur.  Nonlinear  analysis  is  a difficult  topic  and  some  of  the 
points  of  particular  concern  in  finite  element  analysis  of  nonlinear  boundary- 
value  problems  are  related  to  non-uniqueness  and  convergence.  Here  we  wish 
to  examine  the  benefits  of  interweaving  iterative  solution  and  adaptive  mesh 
refinement.  The  finite  element  problem  on  each  grid  now  requires  the  solution 
of  a nonlinear  system  of  algebraic  equations.  Evidently,  it  should  be  more 
efficient  to  refine  the  mesh  as  soon  as  the  solution  iterate  on  a given  mesh 
has  assumed  the  approximate  form  of  a converged  solution,  rather  than  compute 
a fully-converged  solution  on  an  inadequate  grid  prior  to  refinement.  This 
implies  that  one  must  estimate  when  to  cease  iteration  and  refine  the  current 

i 
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grid.  The  relative  and  absolute  changes  of  element  residuals  may  be  used  to 
determine  a stopping  criterion. 

To  explore  this  strategy,  we  consider  a class  of  nonlinear  problems 
arising  in  heat  and  mass  flow  in  chemical  engineering.  The  design  of  effect- 
ive catalytic  reactors,  involving  impregnated  porous  catalyst  in  the  solid 
phase  immersed  in  fluid  reactants,  requires  an  understanding  of  the  transport 
processes  in  the  solid  material.  Consider  the  diffusion  of  heat  and  mass  in 
a catalyst  pellet  where  chemical  reaction  takes  place.  These  problems  are 
characterized  by  nonlinear  reaction  rate  terms,  often  involving  exponentials, 
and  will  exhibit  multiple  solution  states  and  interior-layer  or  boundary- layer 
behavior  for  some  ranges  of  the  reaction  rate  parameters. 

Chemical  engineers  have  developed  high-order  global  approximation  methods 
to  solve  these  "effectiveness  factor"  problems  very  accurately  [5].  However, 
these  global  methods  prove  unsatisfactory  for  the  more  difficult  problems 
where  higher  derivatives  may  be  discontinuous  or  boundary  and  interior  layers 
occur.  To  analyze  these  problems  a C^  orthogonal  collocation  scheme  on  finite 
elements  was  devised  and  applied  to  a representative  "effectiveness  factor" 
problem  with  a boundary-layer  [6].  Problems  such  as  these  are  best  handled 
by  discrete  element  methods,  particularly  if  used  in  conjunction  with  an 
adaptive  refinement  algorithm  that  can  grade  the  mesh  appropriately  into  the 
1 ayers . 

The  solid  phase  consists  of  a porous  base  material  impregnated  with  a 
catalytic  material.  Typically,  the  geometry  consists  of  slabs,  cylindrical 
rods  or  spherical  pellets.  For  these  simple  geometrical  configurations  the 
governing  equation  simplifies  to  a nonlinear  ordinary  differential  equation. 

In  particular,  a first  order,  irreversible,  non-isothermal  reaction  in  the 
catalyst  material  may  be  described  by 


r 
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-JT|-  (x3'1^)'  = f(x,u),  0 < X < 1 (18) 

where  a = 1,2,  or  3 is  the  dimension  and  boundary  conditions  are  u'(0)  = 0 
and,  for  example,  -u'(l)  = B(u(l)  - 1),  parameter  B.  For  f = $ u exp{y(l-l/T)} 
T = 1 + 36  + 3(1  - 6)u(l)  - 3u(x)  and  reaction  rate  parameters  4>  = 14.4, 

3 = 0.02,  y = 20.0,  6 = 50,  B = 250,  the  problem  has  multiple  solutions,  the 

_3 

one  of  interest  possessing  a boundary-layer  profile  of  order  10  near  x = 1. 

As  the  Thiele  modulus  <f>  increases  towards  14.4,  successive  iteration 
on  the  nonlinear  term  becomes  ineffective.  Newton  iteration  is  successful 
in  determining  the  boundary- layer  solution  of  practical  interest.  The  Jaco- 
bian matrix  is  easily  evaluated  as  the  nonlinear  term  contributes  only  to 
the  diagonal  entries  and  to  the  last  column.  The  general  sequence  of  compu- 
tations is  depicted  in  the  flowchart  of  Figure  2,  [2]. 

[Figure  2] 

Applying  the  refinement  algorithm  to  this  problem,  we  again  initiate 
solution  with  a coarse  mesh  of  four  elements.  Results  for  three  refinement- 
iteration  experiments  are  presented  in  Figure  3.  The  parameters  C1 , C2  are 
the  residual  tolerances  in  the  stopping  condition  for  successive  Newton  iter- 
ates. The  graph  demonstrates  that  a single  Newton  iteration  per  mesh  is  most 
efficient  for  this  problem. 

Even  if  the  "final"  graded  mesh  is  employed  directly  with  the  same  start- 
ing guess,  more  operations  (1.64  x io6)  are  required  than  in  the  adaptive 
refinement  solution  (1.25  x 10^).  Moreover,  in  general  we  do  not  have  suf- 
ficient advance  knowledge  of  the  solution  behavior  to  affix  such  a severely 
graded  mesh. 

[Figure  3] 

7.  Potential  Flows 

A class  of  linear  and  nonlinear  flow  problems  that  are  of  considerable 
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practical  importance  in  aerodynamics  arise  in  potential  flow  theory.  We  be- 
gin with  a treatment  of  mesh  refinement  applied  to  incompressible  potential 
flow.  In  particular,  the  local  approximation  near  the  leading  edge  of  a 10 
percent  biconvex  airfoil  is  examined.  The  class  of  problems  is  then  enlarged 
to  admit  compressible  flows,  the  nonlinearity  entering  through  the  density. 
Here  a practical  problem  of  current  interest  is  transonic  airfoil  design. 

The  objective  is  to  devise  techniques  and  efficient  algorithms  for  the  tran- 
sonic flow  regime.  We  select  this  class  of  flows  as  representative  nonlinear 
applications  for  refinement  studies.  The  theory  and  strategies  are  directly 
applicable  to  other  nonlinear  problems  as  indicated  in  the  previous  sections 
and  the  brief  discussions  of  other  flow  examples  in  the  concluding  sections  of 
this  paper. 

8.  Incompressible  Flow 

Conservation  of  mass  together  with  irrotationality  lead  to  the  standard 
potential  equation  for  incompressible  two-dimensional  steady  flows.  A®  = 0 
in  ft  . In  this  instance,  let  us  examine  uniform  flow  about  a 10  percent 
biconvex  airfoil  at  zero  angle  of  attack  • 

The  potential  function  for  the  fully-infinite  exterior  flow  field  may 
be  determined  from  classical  complex  variable  theory  [7], 


$(x,y) 


sinh  b cosh  b 

- X 

2 {(sin  a cosh  b)  + (cos  a sinh  b)“} 


(19) 


where 

a = n (tan  1 " tan_1  S+T> 

b = lln  [{(x+l)2  + y2fy{(x-l)2  + y2}'’*] 

n = | cos'1{(1-62)/(1+6)2} 

a 71 


and 


camber  6 . 
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There  is  a stagnation  point  at  the  leading  edge  and  the  solution  is  less 
' smooth  in  this  neighborhood.  Similar  problems  arise  in  other  applications  in- 

volving corners  or  cracks,  and  local  mesh  refinement  becomes  important.  In 
general,  for  the  Laplacian  on  a domain  ft  with  corner  angle  an  , the  leading 
term  in  the  singularity  of  $ at  the  corner  has  the  form  r 7 sin  e/a  . When 
a > 1 the  corner  is  not  convex  and  the  first  derivatives  of  are  unbounded. 
In  the  present  instance,  since  the  flow  field  is  symmetric  the  corner  angle 
an  is  less  than  n so  that  the  potential  and  velocity  components  (first  de- 
rivatives of  4>  ) are  well  behaved  here.  The  magnitude  of  the  velocity  q 
decreases  to  zero  at  the  stagnation  point  and  then  increases  to  a maximum  at 
the  topmost  point  of  the  airfoil.  In  general,  there  are  fewer  than  1/a  de- 
rivatives at  the  stagnation  point  and  1 + 1/a  derivatives  in  a mean-square 
sense.  For  1/2  _<  a < 1 there  are  only  two  derivatives  in  a mean-square  sense 
so  the  solution  4>(x,y)  is  in  H2.  If  the  airfoil  is  at  a small  angle  of  attack 
the  symnetry  condition  cannot  be  applied  so  that  1 <_ a < 2 and  the  solution  is 
only  in  H^.  In  the  nonlinear  problem  to  follow,  we  are  no  longer  dealing  with 
Hilbert  spaces  and,  moreover,  the  nature  of  the  singularity  and  space  are  not 
known. 

In  the  following  numerical  experiment  we  examine  the  merit  of  refining 
locally  in  a subregion  near  the  stagnation  point.  The  effect  of  this  strategy 
for  nonlinear  analysis  of  compressible  and  transonic  flows  is  then  considered. 

As  the  flow  field  is  symmetric,  it  is  sufficient  to  consider  a single 
quadrant.  The  complex  variable  solution  is  utilized  to  provided  Dirichlet 
data  in  the  far  field  on  a remote  but  finite  boundary  - the  first  "window" 
in  a sense.  Our  mesh  of  triangular  elements  is  generated  automatically  in 
the  region,  as  depicted  in  Figure  4.  Here  the  principles  of  conformal  mapping 
are  introduced:  we  identify  a map  between  a reference  "rectangular"  domain 
and  the  actual  flow  region.  The  mapping  is  defined  by  solving  the  Dirichlet 
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problems  for  the  Laplacian,  V x = 0 and  V y = 0,  on  the  (x,y)  reference  region. 
Finite  element  techniques  are  well-suited  to  this  task  and  the  same  program 
can  be  used  as  in  the  flow  computations.  The  nodal  solution  vectors  of  coor- 

/v 

dinate  values  {x^},^.}  at  nodes  {i}  in  ft  define  the  nodal  locations  in  ft  [8]. 

This  data  basis  and  the  generated  data  set  are  used  repeatedly  with  very 
minor  modifications  for  successively  refined  subregions.  Two  further  refine- 
ments termed  Mesh  2 and  Mesh  3 are  made  in  each  of  the  computations  for  this 
flow  problem.  The  corresponding  subregions  are  marked  in  the  figure. 

[Figure  4] 

8.1  Results 

A finite  element  approximation  $(x,y)  to  the  potential  field  is  calculated 
using  the  given  mesh  of  linear  elements.  Mesh  1.  This  approximation,  the  exact 
solution,  and  error  at  representative  points  on  the  airfoil  are  listed  in  Table 
1,  and  compared  there  with  corresponding  values  obtained  using  Mesh  2. 

Contour  lines  for  the  relative  error  e(x,y)  = {4>(x,y)  - 0(x,y)>/4>(x,y) } 
are  graphed  in  Figure  5.  The  contour  levels  range  from  0.0  in  the  far  field 
to  .0144  at  approximately  the  1/6  chord  position  on  the  airfoil.  As  antici- 
pated, the  contours  are  widely  separated  away  from  the  obstacle  and  dense  in 
the  near  field  where  the  error  function  rises  abruptly  towards  the  airfoil. 

The  subdomain  ft2  for  Mesh  2 is  superimposed  on  the  plot.  The  root-mean- 
square  (RMS)  error  on  ft.|  for  Mesh  1 gives  1 n ( _-j)  = -5.7  for  N = 63  node 
points.  On  the  subregion  ft2  the  error  for  this  mesh  has  ln(F2  = -5.2 
for  N = 25  and  for  (F3  -| ) = -4.8  for  N = 7. 

Subregion  ft2  is  next  considered  and  the  same  mapping  approach  utilized 
to  generate  Mesh  2.  The  boundary  data  is  interpolated  from  the  previous  finite 
element  solution  on  Mesh  1.  As  part  of  the  investigation,  both  exact  and 
prior  finite  element  solutions  were  used  to  provide  boundary  data  on  ft2  . 
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The  results  are  consistent  with  those  predicted  from  the  theory  earlier 
(following  equation  (15)).  The  question  of  error  propagation  from  a boundary 
is  more  interesting  and  difficult  for  the  nonlinear  transonic  problem.  Dis- 
turbances in  the  supersonic  flow  will  propagate  essentially  unchanged  along 
characteristics  whereas  a maximum  principle  holds  here  for  the  linear  problem. 

Solution  and  error  values  at  surface  nodes  for  finite  element  computations 
on  Mesh  2 are  given  in  Table  1.  Again,  the  error  contour  levels  have  the  same 
qualitative  character  as  obtained  using  Mesh  1 on  in  Figure  6.  However, 
they  are  now  less  densely  clustered  and  range  from  0.0  in  the  "new"  far-field 
to  .007  at  approximately  the  1/6  chord  position  on  the  airfoil.  The  error  is 
markedly  reduced  from  that  of  the  previous  calculation:  the  log  of  the  RMS 
error  on  Mes>1  2 is  lntFg  3)  = -6-2  for  N = 63  node  points.  Follow- 

ing the  previous  procedure,  we  identify  the  subregion  and  compute 
ln(F3  2)  ~ "6-°  for  N = 17.  The  actual  choice  of  Oj  is  more  appropriately 
suited  to  the  nonlinear  problem.  As  we  shall  see  in  the  next  section  the 
nonlinear  compressibility  effects  are  most  important  near  the  leading  edge. 

The  final  solution  on  Mesh  3 with  boundary  data  from  Mesh  2 produces 
a similar  error  contour  plot.  The  log  of  the  RMS  error  on  is 
ln(F3  3)  = -6.35  for  N = 63.  The  error  at  the  stagnation  pont  is  .0015.  The 
numerical  experiments  and  results  are  described  in  greater  detail  in  reference 
[9]. 

[Figures  5-6,  Table  1 ] 

9.  Compressible  and  Transonic  Flow 

The  development  of  the  nonlinear  compressible  flow  field  about  an  air- 
foil is  now  examined.  Let  M^  denote  the  uniform  upstream  Mach  number  in 
the  far  field.  An  appropriate  high  Reynolds'  number  Re  Is  assumed  so  that 

the  flow  remains  attached  well  beyond  the  critical  Mach  number  M . For  M 

c 00 

slightly  greater  than  Mc  , local  supersonic  pockets  form  on  the  upper  and 
lower  surfaces  of  the  airfoil.  Supersonic  Mach  numbers  in  the  pockets  are 
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near  unity  and  there  is  no  distinguishable  steady  shock.  As  increases 

further,  a steady  shock  appears  at  the  rear  boundary  of  the  supersonic  pocket 

This  occurs  when  Mmax  is  approximately  1.05.  If  continues  to  increase 

the  supersonic  region  grows  and  the  shock  moves  downstream,  increasing  in 

strength.  A pressure  jump  occurs  across  the  shock,  and  this  eventually  at 

some  M will  lead  to  boundary- layer  separation. 

00 

Finite  element  computations  of  compressible  subsonic  and  slightly  super- 
critical flows  are  next  carried  out  in  conjunction  with  "windowing"  for  mesh 
refinement.  We  repeat  the  sequence  of  computations  on  the  meshes  described 
in  the  previous  section  for  incompressible  aerodynamic  flows.  For  transonic 
flows  the  problem  is  of  mixed  type,  being  elliptic  in  the  subsonic  region 
and  hyperbolic  in  the  supersonic  pocket.  The  mixed  flow  type,  nonlinearity 
and  shock  discontinuities  are  primary  areas  of  difficulty.  We  shall  consider 
only  the  first  two  of  these  in  this  investigation.  In  particular,  accurate 
and  efficient  nonlinear  solutions  are  sought  in  the  vicinity  of  the 
airfoil.  The  exact  solution  is  now  unknown,  and  we  are  unable  to  examine 
errors  explicitly  as  we  did  in  the  previous  linear  problems. 

The  finite  element  formulation  is. first  summarized.  Using  the  Bernoulli 
relation  together  with  the  adiabatic  equation  of  state  in  the  continuity 
equation,we  have  the  transonic  full  potential  equation  [9] 

a * >i  - o (2o> 

where  M is  the  Mach  numoer  of  the  uniform  remote  flow  ana  Cartesian  tensor 
00 

notation  has  seen  employed. 

Equivalently,  we  may  rewrite  the  equation  in  the  form 

U2  - «2)*xx  - 2Vy*xy  ♦ (a2  . «2Hyy  . o 


(21) 


where  a2  = a2  + (U2  - 4>2  - 4>2)  defines  the  local  speed  of  sound,  y is 

the  gas  constant,  and  is  the  uniform  remote  flow.  The  associated  bound- 
ary conditions  in  the  far  field  are 

2 2 

i + Ux  as  x +y  +oo 

and  on  the  airfoil  $>y/<l>x  = fx  where  f(x,y)  describes  the  upper  surface. 

A corresponding  variational  problem  suitable  for  finite  element  compu- 
tations can  be  constructed.  The  functional  is 

1 = J {1  + ^ M2(l  - $2.)}y/y"]  dA  (22) 

A 


Introducing  a linear  finite  element  approximant  into  the  integrand  and 
carrying  out  the  usual  variational  procedures  yields  the  element  contribution 
to  the  nonlinear  finite  element  system 


3f  “ fle  {’  * 


(23) 


where  A is  the  element  area  and  M = L lJ  + L |J  with  $ (x,y)  = . 

C ~ A~  A ■'•y  C “w  ~C 

This  formulation  and  other  approximate  analyses  are  treated  in  greater 
detail  in  reference  [9].  We  consider  mesh  refinement  and  numerical  solution 
of  the  nonlinear  equations  g($)  = 0 obtained  from  equation  (23).  Element 
contributions  to  the  Jacobian  matrix  in  a Newton-Raphson  iterative  procedure 

for  iterate 


can  be  derived  directly  from  (23).  That  is,  J(<^1+1^  - $^)  = -g^ 


i+1,  and  for  each  element 

1 




Me  + (y-D^^H^fM*)!} 
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(24) 


where  8 = (1  + M2(l  - <l>TM<t)i  . 
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In  the  following  numerical  experiments  Newton-iterative  solution  is  com- 
bined with  mesh  refinement  on  successive  windows  adjacent  to  the  leading  edge. 
9.1  Results 

Nonlinear  finite  element  analysis  with  the  nested  meshes  of  Figure  4 is 
now  examined  for  compressible  flow  about  the  biconvex  airfoil.  Solutions 
are  computed  at  a subcritical  incident  Mach  number  M^  = 0.6  and  at  the  slight- 
ly supercritical  Mach  number  M^  = 0.795.  Isomach  contours  are  graphed  in 
Figure  7 for  the  subcritical  flow  (M^  = 0.6)  on  using  Mesh  1. 

[Figure  7] 

In  the  numerical  experiments  we  combine  refinement  of  subregions  with 

Newton  iteration  to  a variable  tolerance  (error  level).  For  example,  the 

_2 

finite  element  solution  on  Mesh  1 is  computed  to  a tolerance  x = 10  in 
the  Newton  iteration.  This  requires  6 Newton  iterates  from  the  incompress- 
ible starting  solution.  The  solution  is  then  sought  on  with  Mesh  2 and 

_3 

to  a tolerance  of  x = 10  in  successive  Newton  iterates.  Finally,  using 

the  solution  from  Mesh  2,  iteration  proceeds  on  subdomain  and  Mesh  3 at 
-4 

tolerance  x = 10  . 

The  behavior  of  the  error  in  successive  Newton  iterates  is  analogous  to 
that  of  the  error  function  in  the  linear  potential  problems  treated  earlier. 

In  Figure  8 contours  of  solution  iterate  error  on  Mesh  1 are  graphed.  Con- 
tour values  range  from  -.006  to  zero.  The  log  of  the  RMS  iterate  error  is 

_3 

-6.4  for  N = 63.  Continuing  to  Mesh  2 on  ^ at  tolerance  10  the  error 

-4  -4 

contours  range  from  -2.6  x 10  to  1.4  * 10  and  ln(e2)  = -14.4  after  2 
iterations.  Values  of  $ and  iterate  error  |e|  at  points  on  the  airfoil 
surface  are  given  in  Table  2.  Finally,  on  fig  with  Mesh  3 and  a starting 
iterate  interpolated  from  the  current  iterate  on  Mesh  2,  at  x = 10"^  the  con- 
tours range  from  -8.1  x 10"7  to  5.4  x 10"7  and  are  shown  in  Figure  9.  The 
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log  of  tne  RMS  iterate  error  is  -14.0  after  3 Newton  iterations. 

The  important  point  is  that  the  absolute  errors  in  each  solution  iterate 
are  near  zero  away  from  the  airfoil  and  increase  to  a maximum  near  the  leading 
edge.  Again,  we  are  exploiting  this  feature  of  the  error  behavior,  but  in 
this  instance  as  applied  to  iteration  errors  rather  than  the  solution  error. 
Rapid  convergence  in  the  far  field  allows  us  to  reduce  the  domain  and  use 
the  previous  "coarse  mesh  solution"  to  interpolate  an  excellent  starting 
guess  on  the  new,  finer  mesh  in  the  near  field. 

We  do  not  mean  to  infer  that  the  actual  error  has  the  same  behavior. 

In  fact,  we  know  that  the  iterate  error  at  the  boundary  is  no  better  than 
that  of  the  previous  mesh.  Yet,  it  may  be  argued  that,  relative  to  this 
error , the  errors  in  the  interior  of  following  subregions  are  being  reduced 
to  zero.  The  final  result  is  a solution  in  the  neighborhood  of  interest, 
the  near  field  of  the  airfoil,  and  on  a fine  mesh  that  is  of  satisfactory 
accuracy.  Furthermore,  this  nonlinear  solution  has  been  obtained  very  ef- 
ficiently. 

[Figures  8-9,  Table  2] 

As  we  remarked  earlier,  the  strategy  proves  useful  if  the  successive 
boundaries  remain  in  the  subsonic  flow.  The  refinement  approach  is  also 
successful  for  supercritical  flows.  At  M^  = .8  there  is  a small  supersonic 
pocket  adjacent  to  the  topmost  part  of  the  airfoil.  On  Mesh  1 we  find  a 
single  element  of  the  coarse  mesh  is  supersonic.  Progressing  to  Mesh  2, 
there  are  several  smaller  elements  that  are  supersonic  and  we  obtain  much 
better  resolution  of  the  sonic  line  and  mixed  flow. 


To  achieve  convergence  of  the  Newton  iteration  at  supercritical  Mach 
numbers,  we  require  an  accurate  starting  iterate  on  an  appropriate  mesh. 
In  the  present  case,  for  mildly  transonic  flows  accurate  computations  of 
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the  flow  field  near  the  leading  edge  can  be  obtained  as  follows:  (1)  Begin 
with  an  initial  mesh  on  the  largest  domain  and  compute  a potential  solution; 

(2)  Using  the  same  mesh  work  up  in  Mach  number  and  with  few  Newton  iterates 
until  the  dominant  changes  in  the  solution  are  localized  near  the  airfoil; 

(3)  Introduce  the  new  mesh,  interpolating  nodal  solutions  from  the  most  re- 
cently calculated  vector  on  the  previous  mesh;  (4)  Continue  iteration  at 
this  Mach  number  on  the  new  subregion  and  on  following  subregions. 

If  accurate  approximations  at  the  higher  Mach  numbers  cannot  be  attained 
on  the  initial  mesh  then  errors  at  newly-introduced  boundaries  may  lead  to 
difficulties.  Alternative  iterative  methods  such  as  point  relaxation  schemes 
may  be  applied  to  make  the  basic  system  solution  more  robust.  These  are  cur- 
rently being  investigated  [8]. 

10.  Other  Problems  and  Refinement  Techniques 

In  this  final  section  a brief  overview  of  other  techniques  and 
applications  is  presented.  Within  this  more  general  context  one  can 
examine  competitive  refinement,  enrichment,  mesh  equidistribution  and 
subregion  refinement  strategies  and  the  classes  of  applications  to  which 
they  are  best  suited.  Directions  of  continuing  research  in  refinement 
theory  and  technique  are  described.  In  particular,  the  implications 
of  these  techniques  to  time-dependent  problems  and  some  related  pro- 
blems peculiar  to  refinement  are  noteworthy. 

10.1.  Refinement  Techniques 

Mesh  refinement  in  one  dimension  is  relatively  straightforward. 


The  techniques  described  earlier  for  mesh  enrichment  guided  by  the 
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residual  are  qualitatively  similar  to  utilizing  the  truncation  error  in 
finite  difference  formulations.  An  alternative  approach,  termed  "equi- 
distribution"  in  the  finite  difference  literature,  cons  ers  a mesh  of 
grid  points  and  redistributes  the  mesh  according  to  an  error  criterion 
such  as  the  residuals  developed  here.  A simple  strategy  is  "doubling 
and  halving"  whereby  the  mesh  size  h is  doubled  in  regions  where  the 
solution  is  better  behaved  and  existing  elements  halved  in  elements 
where  refinement  is  desired.  Other  "monitor"  functions  such  as  arc 
length  and  measured  of  the  local  rate  of  change  of  solution  have  been 
applied,  particularly  in  finite  difference  and  finite  element  colloca- 
tion computations.  A survey  and  bibliography  is  forthcoming  [11]. 

One  may  also  refine  by  increasing  the  order  of  the  element,  a relatively 
easy  process  in  one  dimension  but  obviously  sensitive  to  the  quality  of 
the  initial  mesh. 

Refinement  in  two  or  three  dimensions  is  more  difficult.  Elements 
meet  along  edges  and  surfaces  respectively  and  refinement  will  often 
introduce  new  nodes  on  the  boundary  of  the  original  element.  One  or 
two  internal  refinements  that  do  not  introduce  any  new  nodes  on  the 
old  boundary  might  be  acceptable.  Engineering  experience  has  shown 
that  idealizations  with  slender  elements  may  often  yield  poor  results. 
Mathematically,  we  require  that  the  finite  element  basis  should  be 
uniform.  This  amounts  to  a constraint  on  the  element  geometry.  If  the 
finite  element  approximation  U(x,y)  coincides  with  the  exact  solution 
at  the  nodes  of  a triangular  element,  then 
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where  K is  a bound  on  |3  U/Bx^x^  | , d is  the  diameter  of  the  ele- 
ment and  0 is  the  largest  angle  of  the  triangle.  The  need  for  caution 
is  evident  from  the  dependence  of  the  bound  on  angle  0 . Also,  slender 
triangles  may  lead  to  deterioration  of  the  numerical  condition  of  the 
algebraic  system. 

The  natural  refinement  of  a three-node  triangle  is  to  introduce  new 
nodes  at  the  side  mid-points  and  subdivide  the  triangle  into  four  con- 
gruent subtriangles.  Additional  midside  nodes  are  shared  with  adjacent 
elements.  Across  an  interface  between  refined  and  unrefined  regions 
the  approximation  must  be  constrained  to  ensure  continuity.  The 
constraint  can  be  embedded  directly  in  the  basis,  and  the  approach  has 
been  used  successfully  for  refinement  towards  a singular  point  J12], 
Alternatively,  Lagrange  multipliers  may  be  introduced  to  satisfy  the 
continuity  requirement  on  the  interface  S between  refined  and  unre- 
fined regions.  For  potential  flows  the  variational  problem  becomes: 

make  stationary  the  functional 


I = 


12  2 

+ * )dfl  + 

2 x y 


A(x,y)  (4>+  - 0_)ds 


(26) 


where  and  represent  the  approximations  in  the  refined  and  un- 
refined subregions  adjacent  to  S . The  argument  can  be  applied  using 
penalty  constraints  rather  than  multipliers  in  the  usual  manner  and 
utilised  in  nonlinear  problems  such  as  the  compressible  flow  examples 
considered  earlier.  In  fact,  for  the  transonic  small  disturbance 
equation  [9]  the  modified  functional  is  very  similar  to  that  above. 
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An  approximate  analysis  requires  finite  element  expansions  for  both 
<f(x,y)  in  ft  and  X(x,y)  on  S slightly  increasing  the  size  of 
the  algebraic  problem.  Details  of  the  formulation  and  implementation 
are  described  elsewhere  [13]. 

For  many  applications,  excluding  solutions  possessing  layers  and 
singularities,  only  moderate  refinement  is  required.  In  this  event 
the  additional  complexity  arising  from  the  imposed  constraint  on 
exposed  nodes  may  be  avoided  by  alternative  refinements.  One  scheme, 
that  has  been  utilized  in  an  algorithm  for  linear  and  nonlinar  problems, 
is  to  bisect  the  largest  angle  of  a triangle  and  continue  by  subdividing 
the  adjacent  element  on  the  newly  divided  side  [14].  The  largest  angle 
of  the  second  triangle  may  not  be  bisected  so  that  some  slender  ele- 
ments may  develop.  Yet  for  many  applications  it  is  unlikely  that 
there  would  be  extensive  bisections  of  a small  angle  early  in  the 

refinement  process.  In  [14]  this  approach  is  applied  to  Au  = 1/4  u3/r3 

4 2 

in  a triangular  domain  with  9u/3n  + u - r =0  on  the  slant  side  and 
u * /c  (the  solution)  on  the  remainder  of  the  boundary.  Refinement 
proceeds  from  an  initial  triangulation  of  3 elements  to  60  elements 
after  7 steps.  The  maximum  error  is  3.6  x 10  ^ . 

The  approach  of  refining  nested  subregions  that  is  used  in 
the  potential  flow  experiments  of  earlier  sections  has  been  in  use 
at  least  since  the  mid-sixties. in  finite  difference  computations  and 


probably  was  practised  much  earlier  than  this.  In  industrial  research 
applications  it  has  been  in  use  in  finite  element  computations  since  the 
late  sixties  and,  applications  to  fluid  mechanics  problems  have  been 
made  recently  [15,16].  In  [15]  the  method  is  applied  to  viscous  flow 
through  a two-dimensional  diverging  channel.  The  Reynolds  number  of 
the  flow  is  increased  until  primary  vortices  appear.  In  this  particular 
application  the  separation  region  in  the  channel  is  studied  in  detail  by 
re-meshing  on  the  pertinent  subregion.  Oscillatory  behavior  of  the  solu- 
tion at  moderate  Reynolds  numbers  (around  37.5)  necessitate  the  use  of 
under-relaxation  methods.  A similar  point  arises  in  the  analysis  of 
supercritical  compressible  flows,  but  the  behavior  of  the  iteration 
improves  on  successive  nested  meshes.  Elsewhere  [16]  the  incompressible 
viscous  flow  problem  has  been  studied  numerically  using  a primitive 
variable  formulation  with  nested  refined  subregions.  It  is  observed 
from  experiments  that  the  procedure  again  is  suitable  for  detailed 
examination  of  separation  with  closed  recirculation  regions.  As  in 
the  case  of  the  supersonic  pocket  studies,  details  of  the  flow  near  the 
point  of  interest  are  not  evident  on  the  coarse  mesh  solution  but  are 
resolved  on  subregion  refinement.  For  a rectangular  subregion  the 
velocity  boundary  data  may  be  interpolated  on  three  sides  from  the  pre- 
vious coarse-mesh  solution.  On  the  remaining  side  numerical  experiments 
indicate  that  the  pressure  distribution  is  required  to  ensure  that  mass 
is  conserved  in  the  subregion.  If  velocity  is  also  specified  on  the 
fourth-side  numerical  instabilities  arise  as  the  incompressibility 
condition  is  satisfied  only  in  an  integral  sense  on  the  subregion. 
Detailed  numerical  results  are  obtained  for  corner  recirculation  in  a 
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driven  cavity  flow  at  Reynolds  number  100. 
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Recent  theoretical  work  in  the  form  of  the  residual  analysis  [2,  16] 
is  establishing  a theoretical  foundation  for  existing  and  future  refine- 
ment strategies  in  finite  element  computations  and  there  are  many 
challenging  problems  in  all  areas  of  analysis,  algorithm  and  implementation. 
These  difficulties  notwithstanding.,  in  closing,  we  raise  the  issue  of  time- 
1 dependent  analysis  and  adaptive  refinement.  Equidistribution  techniques 
such  as  the  "doubling  and  halving"  strategy  are  being  applied  to  wave 
propagation  problems.  Consider  a propagating  step-front,  a model 
problem  frequently  studied  of  convection-diffusion  processes.  A fine 
mesh  is  desired  in  the  vicinity  of  the  propagating  front.  A strategy 
currently  being  explored. utilizes  a residual  analysis  within  a semi- 
discrete finite  element  formulation,  enriching  the  spatial  mesh  at  times 
and  locations  determined  by  the  error  criterion. 

11.  Concluding  Remarks 

An  eventual  goal  might  be  to  formulate  and  implement  a theory  for 
adaptive  refinement  whereby  the  analyst  may  specify  conditions  such  as  de- 
sired global  or  local  accuracy,  processing  cost,  or  storage  limitations,  to 
determine  an  appropriate  mesh.  Realistically,  when  one  deals  with  general 
linear  and  nonlinear  problems  this  ideal  can  not  be  attained  - the  degree 
of  difficulty  of  the  problems  and  their  diversity  is  too  severe.  As  we  have 
demonstrated,  substantial  progress  can  be  made  if  one  restricts  the  strate- 
gies to  specific  classes  of  problems  - here  certain  linear  and  nonlinear 
flow  and  transport  problems. 

Even  if  the  more  ambitious  general  goal  is  ignored,  one  can  apply  re- 
finement analysis  and  algorithms  to  great  advantage.  This  is  particularly  true 
if  the  problems  considered  have  layers  or  singularities,  especially  if  their 
location  is  not  well  known  in  advance.  In  treating  nonlinear  problems,  much 
of  the  underlying  computation  required  to  generate  an  accurate  solution 


iterate  and  also  an  appropriate  mesh,  can  be  carried  out  economically  on 
coarse  meshes.  The  nonlinear  iteration  and  adaptive  mesh  refinement  can 
thus  be  interwoven  to  produce  efficient  algorithms. 
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